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Abstract 


Minimal Polynomial Extrapolation (MPE) and Reduced Rank Extrapolation 
(RRE) are two polynomial methods used for accelerating the convergence of se¬ 
quences of vectors {x m }. They are applied successfully in conjunction with fixed- 
point iterative schemes in the solution of large and sparse systems of linear and 
nonlinear equations in different disciplines of science and engineering. Both meth¬ 
ods produce approximations s k to the limit or antilimit of { x m } that are of the 
form Sk = oli x i with X)i=o 7* = h for some scalars 7 The way the two 
methods are derived suggests that they might, somehow, be related to each other; 
this has not been explored so far, however. In this work, we tackle this issue and 
show that the vectors s^ PE and s RRE produced by the two methods are related 
in more than one way, and independently of the way the x rn are generated. One 
of our results states that RRE stagnates, in the sense that s RRE = s RR f , if and 
only if s| IPE does not exist. Another result states that, when s^ PE exists, there 
holds 

= Hk-iSk-i + ^< PE with + v k , 

for some positive scalars fj, kl Hk- 1 , and v k that depend only on s RRE , s RRE , and 
s^ PE , respectively. Our results are valid when MPE and RRE are defined in 
any weighted inner product and the norm induced by it. They also contain as 
special cases the known results pertaining to the connection between the method 
of Arnoldi and the method of generalized minimal residuals, two important Krylov 
subspace methods for solving nonsingular linear systems. 
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1 Introduction 


Minimal Polynomial Extrapolation (MPE) of Cabay and Jackson [5| and Reduced 
Rank Extrapolation (RRE) of Kaniel and Stein [15] , Eddy [6], and Mesina 116] are 
two polynomial methods of convergence acceleration for sequences of vectors]]] They 
have been used successfully in different areas of science and engineering in accelerating 
the convergence of sequences that arise, for example, from application of fixed-point 
iterative schemes to large and sparse linear or nonlinear systems of equations. 

These methods and others were reviewed by Smith, Ford, and Sidi [28] and, more 
recently, by Sidi [24j . Their convergence and stability properties were analyzed in the 
papers by Sidi [20], [23], Sidi and Bridger [25], and Sidi and Shapira [26], [27]. Their 
connection with known Krylov subspace methods for the solution of linear systems of 
equations was explored in Sidi [21] . In Ford and Sidi [9], they were shown to satisfy 
certain interesting recursion relations. Efficient algorithms for their implementation 
that are stable numerically and economical computationally and storagewise were de¬ 
signed in Sidi [22] , Finally, Chapter 4 of the book by Brezinski and Redivo Zaglia [3] 
is devoted completely to vector extrapolation methods and their various properties. 

From the way they are derived, one might suspect that MPE and RRE are somehow 
related. Despite being intriguing and of interest in itself, this subject has not been 
investigated until now, however. In this work, we undertake precisely this investigation 
and show that the two methods are very closely related in more than one way. A partial 
description of the results of this investigation are given in the next paragraph. 

Let {x m } be an arbitrary sequence of vectors in C N endowed with a general (not 
necessarily standard Euclidean) inner product and the norm induced by it, and let 
s fc IPE an d s RRE be the vectors (approximations to lim m _ > , 00 x m , for example) produced 
by MPE and RRE from the k+2 vectors xo, *i, ■ ■ ■, x k+ \. It is known that s RRE always 
exists, but s][f PE may not always exist. One of our results states that, RRE stagnates, 
in the sense that 

s RRE _ s RRE ^ s M p E j oes no £ ex j s £_ (1.1) 

Another result states that, when sJ) IPE exists, there holds 

Vksf RE = Hk-\sf R f + PE with /a k = + v k , (1.2) 

for some positive scalars yL kl fi k -±, and u k that depend only on s RRE , s RRE , and s^ PE , 
respectively. The precise results and the conditions under which they hold will be 
given in the next sections H 

When the sequence {x m } is generated from a linear nonsingular system of equations 
x = Tx+d via the fixed-point iterative scheme x m+ i = Tx m +d , the vectors s^ PE and 
s RRE are precisely those generated by, respectively, the Full Orthogonalization Method 
(FOM) and the method of Generalized Minimal Residuals (GMR)—two important 
Krylov subspace methods for solving linear systems—as these are being applied to 

Whe formulations of RRE given in Kaniel and Stein m and Mesina m are essentially the same, 
but they are entirely different from that in Eddy [6[. The mathematical equivalence of the different 
formulations is shown in Smith, Ford, and Sidi in [281 . 

throughout this work, we will use boldface lowercase letters to denote column vectors. In partic¬ 
ular, we will denote the zero column vector by 0. Similarly, we will use boldface upper case letters to 
denote matrices. 
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the linear system (I — T)x = d, starting with x$ as the initial approximation to the 
solution. This is so provided all four methods are defined using the same weighted 
inner product and the norm induced by itJl 

FOM was developed by Arnoldi [T], who also presented a very elegant algorithm, 
which employs an interesting process called the Arnoldi-Gram-Schmidt process , for 
computing an orthonormal basis for a Krylov subspace. For a discussion of FOM and 
more, see also Saad E3- Different algorithms were given for GMR by Axelsson E],by 
Young and Jea [33], by Eisenstat, Elman, and Schultz j8j, known as Generalized Conju¬ 
gate Residuals (GCR), and by Saad and Schultz [T9] , known as GMRES. GMRES also 
uses the Arnoldi-Gram-Schmidt process, and is known to be the best implementation 
of GMR. For Krylov subspace methods in general, see the books by Greenbaum m, 
Saad HH1, and van der Vorst [29] . 

Now, there are interesting connections between the vectors generated by FOM and 
GMR, and by further Krylov subspace methods, and these connections have been 
explored in Brown [4] and Weiss m originally. This topic has been analyzed further 
in the papers by Gutknecht pf3], pH], Weiss [32], Zhou and Walker [34], Walker [30] . 
and Eiermann and Ernst [7], by using weighted inner products and norms induced by 
them. 

In view of the mathematical equivalence of MPE to FOM and of RRE to GMR when 
{ x m } is generated from linear systems , the results of the present work for MPE and 
RRE [in particular, m and (11.21) ] are precisely those of [4] and [31] in the presence 
of such { x m }. Clearly, our results pertaining to the relation between MPE and RRE 
have a larger scope than those pertaining to FOM and RRE because they apply to 
sequences obtained from nonlinear systems, as well as linear ones, while FOM and RRE 
apply to linear systems only. Actually, our results apply to arbitrary sequences {x m }, 
independently of how these sequences are generated. In this sense, the connection 
between MPE and RRE can be viewed as being of a “universal” nature. We wish to 
emphasize that (i) a priori, it cannot be assumed that MPE and RRE are related when 
applied to vector sequences {x m } arising from nonlinear systems, and (ii) in case there 
is a relationship, it cannot be concluded, a priori, what form it will assume. In view 
of this, the fact that MPE and RRE are related as in m and (|1.2I) in the presence 
of arbitrary sequences {x m }, whether generated linearly or nonlinearly or otherwise, 
is quite surprising. 

The purpose of this work is twofold: 

1. In the next section, we (i) redefine MPE and RRE using a weighted inner product 
and the norm induced by it, and (ii) develop a unified algorithm for their imple¬ 
mentation, thus also providing the theoretical background necessary for the rest 
of this work. 

2. In Section [3] we state and prove our main results showing that MPE and RRE, 
as redefined in Section [2l are closely related. Following this, in Section U we 
discuss the application of our results to sequences {x m } generated from a linear 
nonsingular system of equations via fixed-point iterative schemes, and show that 

3 MPE and RRE were originally defined in C N with the standard Euclidean inner product and 
the norm induced by it. In subsequent work by the author and his co-authors, their definitions were 
generalized by allowing general inner products and norms. The algorithms for implementing MPE and 
RRE given in [22] still use the standard Euclidean inner product and the norm induced by it, however. 
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our main results reduce to the known analogous results of [1] and m that pertain 
to FOM and GMR, also when all four methods are defined using a weighted inner 
product and the norm induced by it. 

The weighted inner product (•, •) and the norm [-J induced by it (both in C^) are 
defined as in 

(y,z) = y*Mz and \z\ = y/(z , z) = Vz*Mz, (1.3) 

where M G C NxN is a hermitian positive definite matrix^ 

For the standard Euclidean inner product and the vector norm induced by it, we 
will use the notation 

( y,z) = y*z and ||z|| = \Jz*z. ( 1 - 4 ) 

A most useful theoretical tool that makes our study of MPE and RRE run smoothly 
is a generalization of the QR factorization of matrices, which we call the weighted QR 
factorization. This version of the QR factorization seems to have been defined and 
studied originally in the papers by Gulliksson and Wedin m and Gulliksson HI], it 
turns out to be the most natural extension of the ordinary QR factorization when 
orthogonality of two vectors y, z £ C N is in the sense (y, z) = 0. For convenience, 
and to make this work self-contained, we discuss it briefly in the appendix. For much 
more, we refer the reader to [12] and HU- 


2 MPE and RRE redefined using a weighted inner 
product 

2.1 General preliminaries 

Let {a? m } be a vector sequence in C N . For the sake of argument, we may assume that 
this sequence results from the fixed-point iterative solution of the linear or nonlinear 
system of equations 

x = f(x), solution s ; x £ and f : C N —> C N , (2.1) 

that is, from 

■®m+l — — 0,1,... , (2-2) 

xq being an initial vector chosen by the user. Normally, N is large and f(x) is a 
sparse vector-valued function. Now, when the sequence {x m } converges, it does so to 
the solution s, that is, linim^oo x m = s. In case {x m } diverges, we call s the antilimit 
of {*m}; vector extrapolation methods in general, and MPE and RRE in particular, 
may produce sequences of approximations that converge to s, the antilimit of {x m }, 
in such a case. 

Let us define the vectors Ui via 

Ui — Xi^-i Xi 1 i — 0,1, ... , (2.3) 

4 Recall that the most general inner product in C N is the weighted inner product that is of the 
form ( y,z) = y*Mz, M being a hermitian positive definite matrix. Of course, in the simplest case, 
M = diagfai,..., a N ) with an > 0 Vi, so that {y, z) = J2? =1 otiyiZi and [zj = oti\zi \ 2 . Finally, 
when M = /, we recover the standard Euclidean inner product and the norm induced by it. 
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and the N x (k + 1) matrices U k via 

U k = [u 0 \ui\ ■ ■ ■ \u k ], k = 0,1,.... (2.4) 

Of course, there is an integer ko < N, such that the matrices U k , k = 0, 1,... ,ko — 1, 
are of full rank, but U ko is not; that is, 

rank (U k ) = k + 1, A; = 0,1,, &o — 1; rank ( U ko ) = ko- (2-5) 

(Of course, this is the same as saying that {ao, tti,..., itfc 0 _i} is a linearly independent 

set, but {uq,ui, ..., u ko } is not.) 

Then, both MPE and RRE produce approximations s k (with k < ko) to the solution 
s of ( 12 . 11 ) that are of the form 

k k 

s k = ^2liXi] 5Z 7i = 1, ( 2 - 6 ) 

2=0 2=0 

for some scalars 7 j. On account of the condition X/f=o 7* = h and because X{ = 
xq + Sj=o u j 1 we can rewrite (12.61) in the form 

k —1 k 

= 7 *’ J = 0,1,.■ • - 1, (2.7) 

j =0 i=i+l 

which can also be expressed in matrix terms as in 

Sfc = x 0 + U k -i£, £ = [£o,£l, • ■ ■ ,^fc-i] T . ( 2 . 8 ) 

We will make use of both representations of s k , namely, (USD and later. 

The 7j and for MPE are, of course, different from those for RRE, in general. 

In the sequel, where confusion may arise, we will denote the vectors s k resulting 
from MPE and RRE by s| IPE and s PR " E , respectively. Similarly, to avoid confusion, we 
will denote the vectors 7 = [70,71, • • • ,7fc] T and £ = [£o,£i, • • • ,£fc-i] T corresponding 
to s k by 7 fc (or 7 ^ PE or 7 PRE ) and £ k (or £^ PE or £ RRE ), respectively, depending on 
the context. 

We now describe how the 7j for MPE and RRE are determined when these methods 
are defined within the context of endowed with a weighted inner product and the 
norm induced by it. 

2.2 Definition of the 7* for MPE and RRE 
2 . 2.1 The 7 i for MPE 

Solve by least squares the linear overdetermined system of equations 


k -1 

^ [ CiUi U k 

i =0 


(2.9) 


for co, ci,..., c k -\. Clearly, this system can be expressed in matrix form as in 

U k -ic' = -u k ; d = [c 0 ,ci,...,c fc -i] T , ( 2 . 10 ) 
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and the least squares problem becomes 


minUL/fc-ic'+ u fe ]. ( 2 . 11 ) 

c! 

Since U k ~1 has full column rank, this problem has a unique solution for c!. Next, let 
Cfc = 1 , and compute 

k 

7 » = —y 1 —j i = 0,l,...,k, provided ( 2 . 12 ) 

Z-ji =0 C * i=0 

From this, we see that s k for MPE exists and is unique if and only if c,; ^ 0. 
(Of course, this means that s k for MPE may fail to exist for some k in some cases.) 

2.2.2 The 7f for RRE 

Solve by least squares the linear overdetermined system of equations 

k 

^ 7 iUi = 0, (2.13) 

i=o 

subject to the constraint 7 ; = 1, for 70 , 71 ,... , 7 k - Clearly, this system too can 

be expressed in matrix form as in 

U k 7 = 0 , 7 = [ 70 , 71 ,... ,7fc] T , (2.14) 

and the constrained least squares problem becomes 

min[[[/fc 7 ]], subject to ££7 = 1; e k = [1,1,... , 1] T G C fc+1 . (2.15) 

7 

(Here e k should not be confused with the kth standard basis vector.) Since U k is of 
full column rank, this problem has a unique solution for 7 . From this, it is clear that 
s k for RRE exists and is unique unconditionally. 

2.3 The special case k = k 0 

With s k for MPE and RRE already defined, we start with a discussion of the case in 
which k = ko. 

Theorem 2.1 Let {x m } be an arbitrary sequence, and let MPE and RRE be as defined 
above. 

1. Provided s^ PE exists, we have s^ PE = s PPE . 

2. Assume the sequence {x m } is generated via (E3D and ( 12 . 21 ) with a linear f(x), 

namely, with f(x) = Tx + d, where T € C NxN is some constant matrix and 
d E C N is some constant vector, and (I — T) is nonsingular. Then s^ PE exists, 
and there holds s^ PE = = s, s being the (unique) solution to x = Tx + d. 

In this case, fco is the degree of the minimal polynomial of T with respect to the 
vector uq . 
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Proof. We start by observing that the matrix U ko - 1 has full rank and that the vector 
u ko is a linear combination of uq, u\, ..., u ko - i- As result, the linear system in (12.101) 
is consistent, hence has a unique solution for c' in the regular sense, this solution being 
also the solution to the minimization problem in (|2.1ip . Letting Ck = 1 and proceeding 
as in (|2.12l) , we obtain the ')'| IPE • A similar argument based on (12.141) and (|2.15l) shows 
that 7 PRE = 7 ^ ipe . This proves part 1 of the theorem. Part 2 can be proved as in 
[28] . for example. ■ 

Since we already know the connection between s]^ PE and s PPE , in the sequel, we 
will consider the cases in which k < ko strictly. 


2.4 Determination of the 7 * via weighted QR factorization 

A numerically stable and computationally economical algorithm for computing the 7 j 
for both MPE and RRE when M = I has been given in Sidi [22]. A nice feature 
of this algorithm is that it proceeds via the QR factorization of the matrices Uk 
and unifies the treatments of MPE and RRE. Of course, in order to accommodate 
the weighted inner product (•, •) and the norm [[•]] induced by it, we need a different 
algorithm. Interestingly, an algorithm that is very similar (in fact, identical in form) 
to the one developed in [22] can be formulated for this case. This can be accomplished 
by proceeding via the weighted QR factorization of Uk- Even though this algorithm, 
just as that in [ 22 ], is designed for computational purposes, it turns out to be very 
useful for the theoretical study of this work concerning the relation between MPE and 
RRE. For some of the details concerning the developments that follow next, we refer 
the reader to [ 22 ]. 

We start with the weighted QR factorization of Uk- For convenience, we describe 
this topic in some detail in the appendix. Since Uk is of full column rank, by Theorem 
IA.1I in the appendix, it has a unique weighted QR factorization given as in 

U k = Q k R k ; Q k € C JVx(fc+1 \ (2.16) 


where Q k is unitary in the sense that Q* k MQ k = I k + 1 since k < N, and R k is upper 
triangular with positive diagonal elements; that is, 


II 

0 * 

Qo 

Ql | 

1 R k = 

r 00 r 0 1 

n 1 

■ ■ r ok 

■ ■ Uk 

? 

(2.17) 






T'kk_ 



q*M qj = 


v i,j] 

r ij = q*Muj 

V i < j; 

ra > 0 

V i. 

(2.18) 


(Note that, having positive diagonal elements and being upper triangular, R k is also 
nonsingular.) Clearly, just as U k has the partitioning U k = [U k - 1 Q k and R k 

have the partitionings 


Q k = [Qk- 1 1 Qk] and R k 


Rk -1 

Pk 

o 1 

T'kk 


Pk = [rok,r lk ,---,r k _ ljk \ T . (2.19) 


We will make use of the following easily verifiable lemma in the sequel: 
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Lemma 2.2 Let 

P e C Nxj and P*MP = Ij. 

Then 

(Py , Pz) = y*z = (y, z) and \Pz\ = \Jz*z = ||z||. 

By this lemma, for arbitrary k, we have 

(QkV,Qk z ) = y* z = (y, z ) and lQk z i = '^ z * z = \\ z W 

and 

l U k zj = \\R k z\\. 

Of these, (12.201) follows from Q* k MQ k = I k +i, while (|2.21l) follows from U k = 
and (12.201) . 


2.4.1 Determination of 7 fc for MPE 

Let us fix c k = 1 and let c = [cq, c\, ..., Cfc] T = 



U fc -i c' + u k = U k c =>- \U k-ic' + u k 1 


Then we have 
= \U k c\ = \\R k c\l 


and the minimization problem in (12.111) becomes, 

min \\R k c\\. 

r' 


By ([2T9D, 


Pk -1 

Pk 


' d ' 


Rk i d + Pk 

0 1 

Tkk 


1 


Tkk 


which, upon taking norms, yields 


||-Rfc c l | 2 — H-Rfe-ic' + p k || 2 + r kk . 

Clearly, by the fact that R k -i is a nonsingular k x k matrix, the minimum of 
with respect to d is achieved when d satisfies 

Rk-ic' + p k = 0 =7 Rk-ic' = -p k =7 d = -R k \p k . 

Of course, d is unique, and so is c. 

With d = — R k ^ 1 p k , the vector 7 fc in MPE is obtained as in 


MPE 

7 k 


c 


e, , c 



Of course, this is valid only when ejc = o°i 7^ O’ hence only when sj! IPE 
The vector c exists uniquely whether s^ IPE exists or not, however. 


( 2 . 20 ) 

( 2 . 21 ) 

Qk R k 


( 2 . 22 ) 

||iifcC|| 

(2.23) 

(2.24) 
exists. 
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2.4.2 Determination of y k for RRE 

By the fact that = ||J2fcj/||, the minimization problem in (12.151) becomes 


min||Rfc 7 ||, subject to ej-y = 1, 

and equivalently, 

nmi'y*(R k R k )-y, subject to ejy = 1. 

By the lemma in |22j Appendix A], the solution for the vector y k in RRE proceeds 
through the following steps: 


R^Rkh = e k , h = [h 0 , hi,..., h k ] T (solve for h). (2.25) 

A = k = tyz (A > 0 always). (2.26) 

z_^i= 0 hi e k h 

yf RE = A h. (2.27) 


Note that h can be determined by solving two (k + l)-dimensional triangular linear 
systems, namely, (i) R* k y = e k for y and (ii) Rkh = y for h. 

For our theoretical study, we need to have "y k in analytical form. This is achieved 
as follows: Substituting h = {R k R k )~ x e k from (|2.25l) in (12.261) and (|2.27l) . we have 


eJ(R* k R k )- ic fc \\R k *e k \\ 2 


RRE _ (R-kR-k) 1 gfc _ R k 1 (R k &k) 

k eJ(R* k R k )-ie k \\R k *e k \\> 
[Here and in the sequel, B~* stands for (R*) _1 = (R^ 1 )*.] 


(2.28) 


(2.29) 


2.5 Unified algorithm for MPE and RRE 

Once the y k have been computed as described above, the computation of s k can be 
achieved via (|2.7p (|2.8I) as follows: First, we compute the vector £ k = [£o>£i> • • • >£fc] T 
via (12.7p . and then, by invoking U k - 1 = Q k _iR k -i, we compute s k via (12.81) . as in 


fc-i 

s k *0 T Q k —i{R k —i£ k ) x 0 + ' rjj 11 j ; 

j =0 

V = Rk-i£ k , y = [r/o,r ? i,...,? 7 fc_i] T . (2.30) 

For convenience, we give a complete description of the unified algorithm in Table 12.11 

















Table 2.1: Unified algorithm for implementing MPE and RRE. 

Step 0. Input: The hermitian positive definite matrix M G C NxN , the integer k, 
and the vectors ®o, *i» * • ■, x k +i- 

Step 1. Compute w* = A Xi = Xi + \ — Xi, i = 0,1,, k. 

Set Uj = [«o | | • • • | Uj\ g C Arx (' 7 + 1 )j, j = 0,1,... . 

Compute the weighted QR factorization of Uk, namely, Uk = Q k Rk', 

Qk = [Qo I Qi I ’ ’ ’ I Qk\ unitary in the sense Q* k MQ k = Ik+i, and 
Rf = [nj]o<i,j<k upper triangular, = q*Muj. 

(Uk -1 = Q k _iRk _i is contained in U k — Qk R k-) 

Step 2 . Computation of ~/ k = [ 7 o, 7 i, • • • , 7 fc] T : 

For MPE: 

Solve the (upper triangular) linear system 

Rf 1 Pkt Pk [^Ok Ulfe , • • • Ufe—l,fc] , C [CO, Cl, • • • , Cfc_i] 

(Note that p k = Q k _ 1 Mu k .) 

Set c k = 1 and compute a = Yli=o c *- 

Set 7 k = c/a; that is, 7 i = Ci/a, i = 0,1,..., k, provided a/0. 

For RRE: 

Solve the linear system 

R* k R k h = e k , h = [h 0 ,hi,..,,h k ] T , e k = [1,1,.. -, 1] T G C fc+1 . 

[This amounts to solving two triangular (lower and upper) systems.] 
Set A = ( Yli =o ^0 • (Note that A is real and positive.) 

Set 7 fc = Ah,; that is, 7 j = A/ij, i = 0,1,..., k. 

Step 3. Compute = [£o,£i,... ,Ck-i] T by 

^o = l-7o; s, s, 1 '■ './• j = 1 ,..., k - 1 . 

Compute s])' IPE and s RRE via 

&k = T Qk —1 (Rk— lCfc) = *0 T Qk—lV- 

[For this, first compute 77 = i 4 -i£fc, h = [ 70 , 7 i> • • •, Vk-iV■ 

Next, set s k = x 0 + 7 Mi-} 
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2.6 Error assessment 


Let us now return to the system of equations in (12.11) . If x is an approximation to the 
solution s of this system, then one good measure of the accuracy of x is (some norm 
of) the residual vector r{x) corresponding to x that is given by 

r (x) = f(x) — x. (2.31) 

This is natural because lim^^g r(x) = r(s) = 0. In case the sequence { x m } is being 
generated as in (12.211 for solving (12.Ill , our measure for the quality of s k will then be 
r(sfe). The following have been shown in (22] : 

• When f(x) is linear [that is, f(x) = Tx + d for some constant matrix T € C NxN 
and constant vector d £ C^], r(s k ) = U k ~t k exactly. 

• When f(x) is nonlinear, U k ~tk serves as an approximation to r(s k ), that is, 
r(sfc) ~ UkJki an d Uk~ik gets closer and closer to r(s k ) as convergence is ap¬ 
proached. 

In addition, J£/fc 7 fc ]], the weighted norm of C/fc 7 fc , can be obtained, without actually 
computing U k J k an d taking its norm, in terms of the quantities provided by the 
algorithm we have just described. This is the subject of the next theorem. 

Theorem 2.3 The vectors C/fc 7 ^ IPE and U fc7 RRR satisfy 

[[C7 fc7fc IPE II = rkkhkl and p7 fc 7 RRE ] = V\. (2.32) 


Remarks. 


1. Of course, ■jk in (12.32(1 is the last component of 7^ IPE corresponding to s^ PE - 
Similarly, A in (12.321) is as defined in (I2.26|) for s RRE . 

2. Clearly, (12.321) is valid for all sequences {x m }, whether these are generated by a 
(linear or nonlinear) fixed-point iterative scheme or otherwise. 


Proof. By (12.211) . we have that [[I7fc7 fc J = \\R k ^ k \\. Therefore, it is enough to look at 
||R fc7 ^ PE || and ||R fc7 RRE ||. 

For MPE, by (|2.22j) . (12.231) . and (|2.24l) . with = l/e p c, we have 


Rk 7f PE = 


e u c 


-( RkC ) = 7 k 


'Y'kk 


— ^kk^k 


Taking norms on both sides, we obtain the result for MPE. 
As for RRE, by (12.291) . we have 


Rk 7 fc RE 


Rf*ek 

R k *e k \ 


2' 


Taking norms on both sides, and invoking (|2.28j) . we obtain the result for RRE. ■ 
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3 MPE and RRE are related 


We begin by restating that since 

U klk = Q k (R klk ) and lU k j k j = \\R k l k \l (3-1) 

and since Q k and R k are the same for both MPE and RRE, the vector that is of 
relevance for both methods is R k l k i an d we turn to the study of this vector. In 
addition, we express everything in terms of the vectors c' and c and the matrices Q k 
and R k , which do not depend either on s^ PE or s RRE . In the developments that 
follow, we will also recall that \\y\\ = \]y y always. 


3.1 R k ~f k for MPE and RRE and an identity 

Assuming that s]^ PE exists, hence ejc / 0, by (12.241) . we first have 


Rnf PE = 


R' k c. 


e k c 


which, upon invoking (12.221) and (12.231) . becomes 

r > „,MPE _ r kk 
Rk^f k „ t 


e u c 


Of course, this immediately implies that 

iii47r E i 

As for RRE, by (12.291) . we have 


Rkif RE = 


Of course, this immediately implies that 

1 


^kk 

I - T i 

Wk c \ 


R k *e k 

\R k *e k \\ 


\\Rklk 


RRE | 


\R k *e k \ 


. _ R ^k 

^ R k G-k — 


RRE 


RRE 112 ‘ 


u k II Rkl k 

We now go on to study R^*e k in more detail. First, by (|2.19p and (12.231) . 

R, 1 = 


r ^-i 

c'/r kk ' 

R~* — 

JD — * 

^k- 1 

0 

Qi 

1 /r kk 

Rk 

. c'*/r kk 

1 /rkk 


Consequently, invoking also e k = 


R k *e k = 


e k -1 


, we have 


r?-* 

R k -1 

0 

_ c!*/r kk 

1 /rkk _ 

r R k 

ll&k-l 


e k -1 


c!*e k _i/r kk + l/r kk 
R k -i^k-i 


G k Gfv kk 


(3.2) 


(3.3) 


(3.4) 


(3.5) 


(3.6) 


(3.7) 
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which, by (13.51) . can also be expressed as in 


1 


Il«fc72 


RRE112 


Rk 7 RRE = 


\Rk-n^ 


RRE 112 




+ 


~ T 

£fc£ 

fkk 


0 

T 


(3.8) 


Clearly, (|3.8D is an identity for RRE relating s?^ E and s PRE ; we will make use of it 
in the developments of the next subsection. (Here t stands for the complex conjugate 
of t.) 

Remark. Recall that vector c exists uniquely for all k < k$. Thus, (|3.8I) is valid 
whether s^ IPE exists or not. 


3.2 Main results 

The following theorem is our first main result, and concerns the case in which sJ) IPE 
does not exist and RRE stagnates. 

Theorem 3.1 1. In case sj) IPE does not exist, there holds 

U k 7^ RE = U k _ l7 RRE . (3.9) 

Consequently, we also have that 


RRE _ RRE 

s k ~ s k -1 • 


(3.10) 


2. Conversely, if (13.101) holds, then sjr PE does not exist. 

Proof. The proof is based on the fact that sj) IPE exists if and only if 0. 

Proof of part 1: Since ejc = 0 when sjf PE does not exist, by (13.81) . 


Rk 7 RRE = 


ll i 47 fe RE || 2 

Taking norms in (13.111) . we obtain 


I*fc-i7£? 


RRE ||2 


RRE 


Rk~nk- 


o 


l|Hrf RE |l = =► Rkt = 


RRE | 


.RRE 


RRE 


r? 

Rk-il k - 


0 


Multiplying both sides of (13.121) on the left by Q k , we have 

U k ~rf RE = Q k (R kl RRE ) = [Q k _i | q k 


td ^.RRE 

Rk-VTk-1 


0 


= Q fc _r(i4_ l7 £5f) = u k . 17^ 


„RRE 


(3.11) 


(3.12) 


(3.13) 


Thus, we have proved the result in (13.91) . 
Rewriting (13.131) in the form 


U kl RRE 


U k 


^,RRE 

'k-1 

o 


(3.14) 
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and invoking the fact that Uk is of full column rank, we obtain 


_ RRE 

7 k = 


-,RRE 

/fc-1 


o 


which, together with (|2.6I) . gives (13.1011 . 

Proof of part 2: By (13. 1011 and (I2.8jl , we have 


RRE i T T z^RRE i j j fRRE RRE 

s k -*0 + ^fc-l€fc - *0 + U k-2$k-l - s k-l > 


from which 


T T £-RRE TT >± 1 ± 1 JD _v T T _ T T 

Uk-lSk — => Uk-l£k ~ U k- 


*.RRE 


RRE 


RRE 


/•HH 

Sfc- 


0 


By the fact that Uk -1 is of full column rank, (13.171) implies that 

‘•RRE i 


£RRE _ 

Sfe — 


Sfc-1 


0 


(3.15) 


(3.16) 


(3.17) 


(3.18) 


which, when combined with the relation [£j = Yli=j+ i 7*> by which, £k-i = 7 k\ hr (12.71) . 
gives (13.151) . Multiplying both sides of (13.151) on the left by Rk, we obtain 


i47fc RE = 




o 


||i4 7 £ RE H = llflfc-17^ 


(3.19) 


Substituting (|3. 1911 in (13.81) . we obtain e E c = 0, and this completes the proof. 


Remark: What Theorem 13.11 is saying is that the stagnation of RRE (in the sense 
that s RRE = s RRE ) and the failure of to exist take place simultaneously. In 

addition, this phenomenon is of a universal nature because it is independent of how 
the sequence {x irL } is generated. 


The next theorem is our second main result, and concerns the general case in which 
*r E exists. 


Theorem 3.2 In case sj) IPE exists, there hold 

1 _ 1 1 
iEopV “ + IvoTT 

and 

UkJf RE Uk- iTgf Uklf PE 

lC/ fc 7^ RE F p7fc-i7j^F Wk 7f PE F' 

Consequently, we also have 

„RRE „RRE „MPE 

^k— 1 ^k 

= I^tEFf + Wl 7 F pe f' 


In addition, 

[[4 7 rre J < 


(3.20) 


(3.21) 


(3.22) 

(3.23) 
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Proof. Since sj) IPE exists, we have ejc 7 ^ 0. Taking the norm of both sides in (13.81) . 
and observing that the two terms on the right-hand side are orthogonal to each other, 
we first obtain 

■ + my, (3 ,4) 


1 


l|i^ RE ll 2 

which, upon invoking (13.31) . gives 

1 


iRk-i^fW 2 


WRkrf 


.RRE112 




+ 


T'kk 


1 


\\Rm% 


.MPE112 


(3.25) 


The result in (13.201) follows from (13.251) and (13.11) . 
Next, invoking (13.21) and (13.31) in (13.81) . we obtain 


1 


\\Rk^ 


RRE 112 


Rk 7 RRE = 




,RRE 112 


td ^.RRE 


+ 


1 


11-^fcTfc 


.MPEII2 


Rk ik 


MPE 


Multiplying both sides of (13.261) on the left by Q k , and invoking (13.11) and 


(3.26) 


Qk 


Rk-nl* 


,RRE 
1 


— [Qk -1 I Qk 


o -,RR 

Rk-llk- 


,RRE 
1 


0 


we obtain (13.211) . 

Let us rewrite (13.211) in the form 


= Qfc-r^fe-iTfe-F) = ^-i7F R F, 

(3.27) 


1 


[E47F re P 


t/^7F RE = 


1 




,RREti2 


_.RRE 
/fc-1 


0 


+ 


1 


C4 7 r E (3-28) 


iu k 7 r E p 


From (13.281) and by the fact that U k is of full column rank, it follows that 


_ RRE _ _ 

\Uk 7 F RE P ” Wk-iT^V 

and this, together with (12.61) . gives (|3.22D . 

Finally, (13.231) follows directly from (13.201) . 

The following facts can be deduced directly from (13.201) : 

lUk 7F RE 1 


_.rre 

Ik-1 


+ 


1 


\Uk 7f PE P 


Ik 


,MPE 


(3.29) 


Wk 7P pe 1 = 


v/i - (tt^7F RE l/[[^-i7F E 


RREJ)2 


when s“ PE exists. (3.30) 


1 

Ic/ fc 7F RE P 


k 


£ 

i&Sk 


1 

wi 7r E p ; 


Sk = {0 < i < k : sP PE exists}. 


(3.31) 


Let us go back to the case in which {x m } is generated as in x m+ \ = f(x m ), 
m = 0 , 1 ,..., from the system x = f(x). As we have already noted, with the residual 
associated with an arbitrary vector x defined as r(x) = f(x) — x , (i) Uk~fk = r ( s k) 
when f{x) is linear, and (ii) U k 7 fc ~ r(s k ) when f(x) is nonlinear and s k is close to 
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the solution s of * = f{x). Then, Theorem 13.21 [especially (13.301) ] implies that the 
convergence behaviors of MPE and RRE are interrelated in the following sense: MPE 
and RRE either converge well simultaneously or perform poorly simultaneously. Let¬ 
ting <f]f PE = [LCfc 7 ^ [PE ]] and 0 RRE = I^fc7 RRE ]], and recalling that </> RRE /(/) RR f < 1 
for all k, we have the following: (i) When </> RRE /<^ RRE is significantly smaller than 1, 
which means that RRE is performing well, (f^ PE is close to ^ RRE , that is, MPE is per¬ 
forming well too, and (ii) when </>^ PE is increasing, that is, MPE is performing poorly, 
</> RRE /0 PPR is approaching 1, that is, RRE is performing poorly too. Thus, when the 
graph of 4>]f PE has a peak for k\ < k < fa, then the graph of </> RRE has a plateau for 
k\ < k < fa. This is known as the peak-plateau phenomenon in the context of Krylov 
subspace methods for linear systems. 


4 Connection with Krylov subspace methods and con¬ 
cluding remarks 

4.1 MPE and RRE on linear systems 

Consider again the linear system of equations x = Tx + d, where the matrix (I — T) 
is nonsingular, and generate {x m } via x m+ i = Tx m + d, m = 0,1,..., with some 
initial vector xq. Apply MPE and RRE to {x m } to obtain the vectors s as before. 
As already stated, Uk~fk = r k = r(sjfc), where r(x) = (Tx + d) — x is the residual 
vector for the system (/ — T)x = d associated with x. In this case, we have the next 
theorem as a corollary of Theorems 13.11 and 13.21 


Theorem 4.1 Let the sequence {x m } be generated recursively via x m+ i = Tx m + d, 
m = 0,1,..., the matrix (I — T) being nonsingular. Let also r(x) = Tx + d — x be the 
residual vector corresponding to x. Let k o be the degree of the minimal polynomial of 
T with respect to uq = x\ — xq . Then, for k < k o, the vectors s^ PE and s RRE obtained 
by applying MPE and RRE to {* m } and their residual vectors r(s^ PE ) = r^ [PE and 
r ( s RR E ) _ r RRE sa f/ LS j‘y following for this special case: 


-j C.RRE _ 


if and only if sf 


MPE 


fails to exist. 


2. In case sP PE exists, there hold 


1 


•.RRET12 


1 


„RREJ2 


+ 


1 


„MPE]12 ' 


r 


RRE 


„RREti2 


Consequently, we also have 


„RRE 
k -1 


K*_i E P 


„MPE 


+ 


pr E F 


-.RRE 


[r RR E J2 


-.RRE 

5 /c-l 


I7- R .T 


,MPE 


+ 


-MPET12 ' 
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In addition, 


„RREti 


< \\r 


RREti 
k -1 1 


(4.1) 


(4.2) 


(4.3) 

(4.4) 
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3. sf™ = s RRE = s, where s is the solution to (I — T)x = d. 


In view of (14. 1 j) . the results in (13.301) and (13.311) become 


and 





I r R RE ) 


(Ir R RE ,/ [r RREj)2 


when sf PE exists 


1 


k i^ PE i 2; 


Sk = {0 < i < k : *r E exists}. 


(4.5) 


(4.6) 


4.2 Equivalence of redefined MPE and RRE to Krylov subspace 
methods for linear systems 

Theorem 2.4 in m concerns the mathematical equivalence of vector extrapolation 
methods to Krylov subspace methods, when all these methods are defined using the 
standard inner product (•, •) and the standard norm || • || induced by (•, •): This the¬ 
orem states specifically that MPE and RRE are equivalent to, respectively, the full 
orthogonalization method (FOM) of Arnoldi and the method of generalized minimal 
residuals (GMR) when 

• MPE and RRE are being applied to the sequence {x m } obtained via x m+ i = 
Tx m + d, m = 0,1,... , with some xq, and 

• FOM and GMR are being applied to (/ — T)x = d, starting with the same initial 
vector xq. 

As stated in Theorem 14.21 below, this theorem holds true also when MPE, RRE, FOM, 
and GMR are defined using the weighted inner product (•, •) and the weighted norm 
H induced by (•,•). 

For the solution of a nonsingular linear system Ax = b, whose solution we denote 
by s, FOM and GMR construct their approximations Wk as follows: Define the residual 
vector corresponding to x by r(x) = b — Ax and denote tq = t{xq) for some initial 
vector xq. Let JCk(A;ro) = spanjro, Atq, ..., A k ~ 1 ro}. Then, for each method, the 
approximation Wk to s is of the form Wk = xq + y such that y £ /C^(A; ro), and y is 
the vector to be determined. Using the weighted inner product (•, •) and the norm [[-J 
induced by it, these methods can be redefined as follows: 

• For FOM, y is determined by requiring that {z, r POM ) = 0 for all z £ ICk{A; r o), 
where r*f 0M = r(m POM ). 

• For GMR, y is determined by requiring that = m.\n. y& ^ k ^. ro ^ Jr (* 0 +?/)]], 

where = r(m)? MR ). 

Then we have the following generalization of Theorem 2.4 in [21] : 
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Theorem 4.2 Consider the nonsingular linear system (I — T)x = d. Apply FOM 
and GMR to this system starting with some initial vector xq. Apply MPE and RRE 
to the sequence {x m } obtained from x m+ i = Tx m + d, m = 0,1,, with the same 
initial vector xq. Then 

W F° M = S MPE and ^GMR = S RRE^ ( 4 .7) 

when all four methods are defined using the same weighted inner product {■, •} and the 
norm [•]] induced by it. Consequently, all of the results of Theorem \4.1\ apply verbatim 
to the vectors -u; EOM and 

Proof. The same as that of [211 Theorem 2.4]. ■ 

These results for FOM and GMR are not new. They were given originally by Weiss 
m and by Brown []4J, and developed further in the papers mentioned in Section [1] 

Note that the vectors u; EOM and -ic)f MR can be obtained numerically by modifying 
the known algorithms for FOM and GMR via the introduction of the weighted QR 
factorization in the Arnoldi-Gram-Schmidt process, in exactly the same way we have 
approached the problem of computing s^ PE and s RRE , the rest being simple. 

Appendix: Weighted QR factorization 

A.l Gram—Schmidt process and existence and uniqueness of the 

weighted QR factorization 

Gram-Schmidt orthogonalization, also called the Gram-Schmidt process (GS), is a pro¬ 
cedure that takes an arbitrary set of linearly independent vectors {ai,...,a s } and 
constructs a set of vectors {q 1 ,..., q s } out of them, such that the q t are orthonormal 
with respect to some specified inner product. In addition, the q i are generated in a 
way that ensures 


span{oi, ...,aj} = span{q 1; ..., qfi, j = 1,..., s. 
This, of course, implies 


a i = r u q 1 

a 2 = r 12 qi + r 22 q 2 

«3 = n 3 q 1 + r 23 q 2 + r 33 q 3 

a s = r 3s q l + r 2s q 2 4 - \~r ss q s , (A.l) 

for some suitable scalars r^. 

Remark: We observe that, with ai,q i E C m , the relationships in (1A.1I) . with no 
conditions imposed on the qj and rij, can be expressed in matrix form as in 

A = QR, (A.2) 
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where 

and 


A = [a\ | d 2 | • • • | a s 1 € C r ‘ 


(A.3) 


Q = [gil9 2 l ••• l</Jer xs , = 


rn H2 • • • r-is 

^22 • • • r 2s 


G C sxs . 


(A.4) 


By this, we mean that (1 A.2I) (I A.4|) follow directly from (IA.1I) . in which the vectors 
q ; and the scalars rjj are not required to satisfy any specific conditions except (IA. Ill 
itself. 

Let (•, •) be a weighted inner product in C m , given as 

(y,z) = y*Gz, G € C mxm hermitian positive definite, 
and let [•]] be the norm induced by (•, •), given as 


M = sj (z,z) = Vz*Gz. 


We now return to the determination via GS of the q ; and the rjj such that the Qj 
satisfy the orthogonality condition ((/,-, qj) = dij. The first of the relations in (1 A. 1 1) 
gives rn and q x as follows: 


i"u = Jai], Qi = a-i/rn. 

With q 1 determined, and by invoking the requirement that ( q\,q 2 ) = 0, the second 
relation is used to determine r\ 2 , r 22 , and q 2 as in 

r 12 = (Qi, a 2 ), a 2 = a 2 - n 2 q x , r 22 = [a 2 1, q 2 = a 2 /r 22 . 

We now proceed in the same way and determine q 3 ,...,q s and the rest of the For 
example, with q 1 ,..., qj_ l already determined and (q { , q,^) = 5 a* for 1 < i, i' < j — 1, 
and by invoking the requirement that ( qi,qj) = 0 for i < j, from the relation dj = 
YH=i r ij ( li i n 11 A. 1 1) . we determine rij for i < j and as follows: 

3 - 1 

rij = (q t , cij ), 1 < i < j; cij = a ? - ^ rjj = [[ajJ, q :j = dj/rjj. 

i =1 

Note that, since the a,; are linearly independent, the vectors dj are nonzero and hence 
rjj > 0 for all j. Note also that r^q, is the projection of dj along q^, i < j. 

Since the columns of the matrix Q are orthogonal with respect to a general inner 
product in C m , and since the most general inner products are necessarily weighted 
inner products, we will call the QR factorization of A given in (|A.2j) (IA.4[) a weighted 
QR factorization. The next theorem summarizes the existence and uniqueness of the 
weighted QR factorization. 
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Theorem A.l Let 


A = [ ai | a 2 | • • • | a s ] G C mxs , m > s, rank( A) = s. 


Let also G G C mxm be hermitian positive definite and define the weighted inner product 
(■, ■) via (y, z) = y*Gz. Then there exist a matrix Q G C mxs , unitary in the sense 
that Q*GQ = I s , and an upper triangular matrix R G C sxs with positive diagonal 
elements, such that 

A = QR. 


Specifically, 


Q 


[Qi\Q2\ ••• \Qs], 



r 12 
T22 


r\ s 

r2s 

r ss 


(Qi, Qj) = Q*Gq, = Sij V i,j, 
r ij = (9i) a j) = QiGdj V i < j; > 0 Vi. 


In addition, the matrices Q and R are unique. 


Proof. Clearly, the columns of Q and the entries of R are precisely those in (1A.1I) . 
as provided by the Gram-Schmidt orthogonalization we have just described, with the 
inner product (•,•). We have thus shown the existence of the QR factorization for A. 
[Recall also the remark following (1A.11) .] 

We now turn to the proof of uniqueness. Suppose that A = Q 1 R \ and A = Q2R2 
are two QR factorizations. Then, by Q\R\ = Q2R2 , we have both R\R 2 1 = Q\GQ 2 
and RiRi 1 = Q\GQi- Since R\R 2 1 and R2R1 1 are both upper triangular with 
positive diagonals, so are Q\GQ 2 and Q^GQ 1 . But Q\GQ 2 = {Q^GQfi}* , and hence 
Q\GQ 2 is lower triangular as well; therefore, Q\GQ 2 is diagonal with positive diagonal 
elements. So is Q* 2 GQ 1 . Now, Q\GQ 2 is unitary since 

{Q\GQ 2 )*{Q\GQ 2 ) = (Q* 2 GQ 1 )(QIGQ 2 ) = (R 2 R^ 1 )(RiR 2 1 ) = Is. 

Since the only diagonal unitary matrix with positive diagonal elements is I s , we have 
Q\GQ 2 = I s - Similarly, we have Q^GQ 1 = I s . Therefore, 

R\R 2 X = I s = RiRi 1 , hence Ri = R 2 . 

As a result, 

Qi = AR 1 1 = AR 2 1 = Q 2 . 

This completes the proof. ■ 

Remark. The weighted QR factorization of Theorem IA.1I reduces to the standard 
QR factorization when G = I m , that is, when (•, •) is the standard Euclidean inner 
product (y, z ) = y*z. In this case, Q is a unitary matrix in the regular sense, that is, 

Q*Q = I S . 
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A.2 Weighted QR factorization by the modified Gram-Schmidt pro¬ 
cess 

It is well known that, in floating-point arithmetic, the orthogonalization process [with 
the standard inner product (y, z ) = y*z\ via GS is very inaccurate, and that this 
problem can be fixed to a large extent by using the modified Gram-Schmidt process 
(MGS) in order to overcome numerical instabilities. In view of this, it seems reasonable 
to suggest that MGS would be preferable to GS with the weighted inner product 
( y,z ) = y*Gz too, even though the (numerical) end result may be influenced by the 
matrix G. Here are the steps of MGS: 

MGS Algorithm 

MGS1. Set rn = [aij and q 1 = ai/rn. 

MGS2. For j = 2,. .., s do 

Set = dj. 

For i = 1,..., j — 1 do 

Compute r t j = (q i ,ap) and compute a y +1 ' = ap — r^q^ 

end do (i) 

Compute rjj = |a^]] and set q t = ap/rjj. 

end do (j) 

It is easy to show that, in exact arithmetic, the scalars and the vectors qj 
computed by MGS are the same as those computed by GS. Comparing the GS and 
MGS algorithms, we also see that their costs are identical. However, they differ in the 
way they compute the [namely, r t j = (q t . afi for GS, while = ( q i , ap) for MGS, 
the rest of the computations being identical], and this may have a hopefully beneficial 
effect on the rest of the computation. For example, once r^q 1 , the projection of a 3 

along q 1 , has been computed, it is immediately subtracted from ap = 03 to give a!p 

( 2 ) 

and then r 23 is computed via (q 2 , a 3 ) instead of (q 2 , 03 ), even though they are equal 
mathematically speaking. 

Before closing, we mention that weighted QR factorizations are discussed in great 
detail by Gulliksson and Wedin m and Gulliksson PH- These papers discuss the use 
of MGS and also of suitable generalizations of Householder reflections in computing 
the weighted QR factorization. Numerical aspects of this computation for diagonal 
matrices G = diag(ai,..., a s ), a* > 0 V i, are studied in [H] for different relative sizes 
of the a*. 
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